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Abstract 

A recently proposed method for computer simulations in the isothermal- 
isobaric (NPT) ensemble, based on Langevin-type equations of motion for 
the particle coordinates and the "piston" degree of freedom, is re-derived by 
straightforward application of the standard Kramers-Moyal formalism. An 
integration scheme is developed which reduces to a time-reversible symplectic 
integrator in the limit of vanishing friction. This algorithm is hence expected 
to be quite stable for small friction, allowing for a large time step. We discuss 
the optimal choice of parameters, and present some numerical test results. 
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I. INTRODUCTION 



Molecular Dynamics (MD) simulations [|]-|3| are a very efficient tool to study the statisti- 
cal properties of thermodynamic systems, especially at high densities where the acceptance 
rates of standard Monte Carlo (MC) simulations [|J are small. They are also very well suited 
to study dynamical properties. 

MD simulations are most naturally performed within the micro canonical (NVE) ensem- 
ble, while the simple Metropolis MC algorithm leads to the canonical (NVT) ensemble. 
However, it is often desirable to study a system within a different ensemble, and hence 
methods have been developed to extend both MC and MD methods to practically every 
ensemble of thermodynamics 0-0- For MC methods, this task is relatively straightforward: 
One starts from the partition function or equilibrium probability distribution function for 
the pertinent degrees of freedom, and defines a Markov process in the space of the latter. 
After verifying the condition of detailed balance, and making use of ergodicity (which often 
can be safely assumed, and in some cases even be shown rigorously), one has established 
that the process will ultimately produce fluctuations within the proper equilibrium. For 
example, for a simulation in the isothermal-isobaric (NPT) ensemble, the system is studied 
in a box with periodic boundary conditions, whose size is allowed to fluctuate. In order 
to keep the system homogeneous, the particle coordinates are instantaneously adjusted to 
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these box fluctuations, via simple rescaling. We shall here be concerned with the simplest 
case only, where the box is just a cube of size L in each direction. By writing 



fi = Lsi (1) 

one introduces reduced coordinates s*j in the unit cube instead of the original coordinates fi 
of the N particles in the box of volume V = L d , d denoting the spatial dimensionality. The 
abovementioned adjustment of the particle configuration to the box fluctuations is facilitated 
by updating L and the independently. The partition function appropriate to the NPT 
ensemble is then 

/■oo r- p 

Z = J dV J d d n... J d d r N exp (-(3U - (3PV) 

POO P P 

= J dV J d%...J d d s N V N exp (-/3U - (3PV) , (2) 

where U denotes the system's potential energy, P is the external pressure, and (3 = 1/(&bT). 
From this one immediately reads off that one has to run a standard Metropolis algorithm 
on the state space (L, {<?«}), using an effective Hamiltonian 

U eff = U + PV- Nk B T In V. (3) 

The MD approach mil to non-microcanonical ensembles [HH^,§-0|, pioneered by An- 



dersen ||, Nose H and Hoover ||10|| , is slightly more involved. Like in MC, one defines 
an additional dynamical variable whose fluctuations allow to keep the thermodynamically 
conjugate variable fixed. In our example, this variable is L, while P is fixed. However, 
the dynamics is not specified via a random process, but rather by Hamiltonian equations 
of motion in the extended space. This requires the definition of canonically conjugate mo- 
menta of the additional variables, and in turn the introduction of artificial masses, which 
have no direct physical meaning but are adjusted in order to set the time scale of the new 
fluctuations. The analysis then proceeds by assuming ergodicity in the extended space, such 
that the statistics is described by a microcanonical ensemble in that space. If the equations 
of motion have been chosen properly, then the equilibrium distribution function of the al- 
gorithm, which is obtained by integrating out the artificial momenta, must coincide with 
that of the desired ensemble. For constant-pressure simulations, often the picture of a "pis- 
ton" is employed. In the last few years many contributions on extending constant-pressure 



methods have been made, treating cases of isotropic boxes as well as non-cubic cells fTT|-|16 
We will here however be only concerned with the case of an isotropic box, as in the original 
Andersen method |J which produces an isobaric-isenthalpic (NPH) ensemble. 

Stochastic dynamics (SD) in its classical form 0] can be viewed as a simulation method 
which is somewhere between MC and MD, sharing with the former the stochastic element 
(which also ensures the ergodicity of the method) and the generation of a canonical (NVT) 
ensemble, while being based on continuous equations of motion, and momenta, like the 
latter. Instead of Hamiltonian equations of motion one solves Langevin equations and adjusts 
the temperature via the balance between friction force and stochastic force (fluctuation- 
dissipation theorem). It is not surprising that this approach can be combined with the 
Andersen method in order to produce an NPT ensemble. It is the purpose of this paper 
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to show that the algorithm can be derived in a very simple and straightforward manner, 
by exploiting the well-known equivalence of the Langevin equation with the Fokker-Planck 
equation |Tj|, and avoiding the complicated reasoning of a recent publication The 
derivation is outlined in Sees. [0], [H| and [PV| , where we treat a straightforward generalization 
of SD to arbitrary Hamiltonian systems, and apply this to the original Andersen NPH 
method, which is thus modified to an NPT method. We then turn to the question of 
numerical implementation. Since SD reduces to standard Hamiltonian dynamics in the 
limit of zero friction, and practical applications are often run for rather small friction, it is 
useful to use an algorithm which reduces to a time-reversible symplectic integrator (TRSI) 
in the zero-friction limit. It is well-known that TRSIs are particularly well-suited for 
Hamiltonian systems PJ^, |l5| . |r^ .^0|], since, except for roundoff errors, they do not prefer 
a particular direction of time and are hence very stable. Our algorithm is derived and 
tested in Sec. [V|, where we apply the method to a simple model system of Lennard- Jones 
particles, and pay particular attention to the question how the parameters should be chosen 



for optimum performance. Finally, we conclude in Sec. VI 



II. GENERALIZED STOCHASTIC DYNAMICS 

Our starting point is a set of generalized coordinates qi and canonically conjugate mo- 
menta Pi, such that the Hamiltonian equations of motion read 

an . m 

dpi dqi 

where Ti is the Hamilton function. These equations of motion are generalized to their 
stochastic versions by adding friction forces with damping parameters 7« ({<&}) and stochastic 
forces with noise strengths cr^ {{qi}) (note that a dependence on the generalized coordinates 
is permitted): 

. m .mm 

q *=dp- ft = -^-*flk ( 5 ) 

where the Gaussian white noise rji satisfies the usual relations 

(tfc(i)) = (Vi{t)Vj{t')) = 2M (t - • (6) 

Equivalently, the stochastic process is described by the Fokker-Planck equation |18[ 



which describes the time evolution of the probability distribution function $ in the full 
space of stochastic variables. For an arbitrary set of variables x v (in our case both q^ and 
Pi) it reads 

where the right hand side of the equation defines the Fokker-Planck operator Lpp. Drift 
and diffusion coefficients and D&J are related to the short-time behavior of the process 
via the Kramers-Moyal expansion [|T^] : 
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D$ = Km ^(Ax,(At)Ax v (At)). (8) 

Straightforward evaluation of these moments for the present case yields directly the Fokker- 
Planck operator: 

i dpi dqi Y % ^ i ^ V ^ a% dpi) ' ' 

A canonical ensemble is generated if the Boltzmann distribution is the stationary dis- 
tribution of the process (due to the stochastic element, the process is usually ergodic, such 
that only one stationary distribution exists). For the present case one finds 

zL FP exp {-pn) = £ A ( 7 . _ pa*) |^exp (-(3H) , (10) 
which vanishes if 

a* = fc fl T 7i . (11) 

This relation is the generalized fluctuation-dissipation theorem which friction and noise have 
to satisfy in order to generate a canonical ensemble. 



III. ANDERSEN EXTENDED SYSTEM 

The Andersen method || uses the box volume V = L d and the scaled coordinates Sj, see 
Eqn. |l], as degrees of freedom. For the particle velocities one obtains 

ri = Lsi + Lsi, (12) 

however, the second term is deliberately omitted in order to achieve independent fluctuations 
of L and Sj. Hence the method amounts to postulating the Lagrangian 

£ = E |W* 2 - E % (A W) + f^ 2 - p ^ (is) 

where denotes the mass of the zth particle, Q is the artificial piston mass or box mass, 
and V{j is the interaction potential between particles % and j (the generalization to three- 
and many-body forces is straightforward). The Hamiltonian is then obtained via Legendre 
transformation 
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where we have used the canonically conjugate momenta 

dC ■ _ dC 



Tl v = ^ = QV, n i = ^9r = mL 2 Si, (15) 

oV dsi 



such that the Hamiltonian equations of motion read 

• 1 



L 2 m, 



■7Ti 



7?i = Lfi (16) 



7 = ^Tl v tl v = V-P, 

where fa is the force acting on the ith particle, and V abbreviates the "instantaneous" 
pressure 

/y being the force acting between particle i and j, while % = Si — Sj. 

Obviously, 7i is a constant of motion. Apart from the kinetic energy of the piston, and 
the deviation of the simulated molecular kinetic energy from the true kinetic energy, which 
are both small corrections [^], 7i is just the enthalpy. For this reason, the method produces 
the NPH ensemble. 



IV. LANGEVIN EQUATION FOR THE NPT ENSEMBLE 



The idea of Feller et al. \TB] was to replace the canonical equations of motion (Eqn. |TB|) 
by a Langevin stochastic process. It was designed to avoid oscillations of the box volume 
("ringing" of the box). In Ref. |l]| an infinite set of harmonic oscillators coupled to the 
box piston was used to prove the correctness of the approach. However, since the original 
Andersen method is based on a Hamiltonian system, the results which have been derived in 
Sec. |II| can be directly applied. Therefore the stochastic equations of motion read 



7?i = Lfi - -jY—Ki + J 'k B T^ifji{t) (18) 

±j 171 j 



ti v = v~p- ^n y + ^k B T lv fj v {t) 
(v?) = (w) = o 

( Vv (t) Vv (t')) = 25(t-t') 

fo?(*V(0> = °' 

while the equations of motion for and V remain unchanged. Here a and (3 denote Cartesian 
directions, and the fluctuation-dissipation relation has already been taken into account. 

This method generates the correct NPT ensemble, as is seem from writing down the 
partition function which naturally arises from the algorithm (cf. Sec. 0): 
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Z = J <m v J d d it x ...J d d 7t N JdV Jd d S!...J d d s N exp (-f3H) , (19) 

where one obviously has to use the Hamiltonian of the Andersen method, see Eqn. [14]. Inte- 
grating out the momenta, one sees directly that this is, apart from unimportant prefactors, 
identical to Eqn. 0, i. e. the correct partition function of the NPT ensemble. 

Note that there is still considerable freedom in the choice of the friction parameters. 
Feller et al. chose ji = and = const.. This is tantamount to coupling only the piston 
degree of freedom to the heat bath. Since this degree of freedom is tightly coupled to all 
the others, the method produces the same NPT ensemble as the more general case 7j ^ 0. 
However, we view this latter case, which also includes a direct coupling of the particles to the 
heat bath, as more advantageous, not for fundamental reasons, but rather for practical ones: 
As in standard SD, every degree of freedom is thermostatted individually. For this reason, 
local instabilities, arising from discretization errors, are efficiently corrected for, without 
spreading throughout the system. Loosely spoken, a particle which happens to be too "hot" 
will be "cooled down" by its local friction, while in the opposite case it will be "heated up" 
by the noise. For this reason, the SD version with 7« 7^ allows for a slightly larger time 
step than the pure Andersen method, while the Feller et al. version (7; = 0, 7y 7^ 0) is not 
more stable, involving only global thermostatting. We have chosen 

li = 7o£ 2 , (20) 
while choosing a constant value for 7y. Then the Langevin equation for 7fj is rewritten as 



7?< = L U ~~ + ^k B T lo mj ; (21) 

this fits naturally to standard SD, where the friction force is —■yoPi/m i = — ^Q-Ki/ (Lrrii). 
Further details on the implementation will be given in Sec. [V]; there we will also discuss the 
choice of the parameters 70 and 7^, as well as of the piston mass Q, which are all irrelevant 
for the statistical properties of the system, but of great importance for the equilibration 
properties of the algorithm. 



V. IMPLEMENTATION 
A. Symplectic Integrator 

Symplectic time-reversible integrators are known to be extremely useful for molecular 
dynamics simulations of Hamiltonian systems P,[7| JT5|JT9| -f22| . This is so because, except for 
roundoff errors, they conserve the phase-space volume exactly (note the intimate relation to 
entropy production), and do not mark any particular direction of time (note that a global 
drift in the algorithm breaks this time-reversal symmetry). Hence they are very stable and 
allow for a large time step. The most common example (actually, the lowest-order scheme) 
is the well-known Verlet algorithm in its various formulations [l|]. 

A particularly transparent way to construct these algorithms |19[ is based on the Liouville 



operator iL, which is just the Fokker-Planck operator in the special case of vanishing friction 
(see Eqn. ||). Noticing that the formal solution of Eqn. ^ is just 
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$(t) = exp(iLt)$(0), (22) 

one focuses directly on the operator exp (iLt) , and decomposes L into a sum of simpler 
operators, 

n 

iL = J2*L k , (23) 
fc=i 

where the imaginary unit is extracted for convenience; L as well as each of the is self- 
adjoint. Now, the exact time development within a time step At is replaced by an approxi- 
mate one by using the factorization 

e iLAt _^ e iL 1 At/2 e iL 2 At/2 e iL n -i At/2 e iL„ At^L^ At/2 ^L 2 At/2 ^LiAt/2 ^4) 

This scheme is automatically phase-space conserving, since each of the operators is unitary, 
and time-reversible symmetric, since the inverse operator is just the original operator, eval- 
uated for —At. Now, if each of the operators Lf. is simple enough, the action of e lLkt on 
a phase space point can be calculated trivially, such that the algorithm is a succession of 
simple updating steps. 

Specifically, for the Andersen equations of motion (see Sec. |TT1|) we choose the operators 

i£, = -£«|- (25) 
iL 2 = -{V-P) ° 



iL 



Rv_d_ 
Q dV 

resulting in the following updating scheme: 

1. j?i(t) -> 7T«(* + At/2) = 7Tj(t) + L(t)fi(t)At/2 

2. Tl v (t) -> U v (t + At/2) = Tl v (t) + (V-P) At/2 

(note that for the evaluation of V, one has to take the old positions s$(t) and the old 
box size L(t), but already the updated momenta 7fj(£ + At/2)) 

3. V(t) -> V(t + At/2) = V(t) + Q- l Il v {t + At/2)At/2 

4. - s t (t + At) = m + Jl^l At 

5. 1/(t + At/2) -> V(* + At) = V(t + At/2) + Q- 1 ^^ + At/2) At/2 

6. rb/(t + At/2) -> U v (t + At) = n y (t + At/2) + (V-P) At/2 

(for evaluation of V, one uses s*j(t + At), L(t + At), and jti(t + At/2)) 
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7. 7T;(t + At/2) -> 7?i(t + At) = 7?i(t + At/2) + L(t + At)/i(t + At) At/ 2 . 

It is often convenient to formulate the algorithm in terms of the conventional variables 

fi(t) = L(t)si(t) pi(t) = Ht)- 1 ^), (26) 

which are however not canonically conjugate with respect to each other. The pressure, in 
terms of these variables, is written as 

and the updating scheme, which now involves various rescaling steps, proceeds as follows: 

1. # = #(t) + /7(t)At/2 

2. V is evaluated using Eqn. [27] with r*j(t), L(t), and pj; then ITy is updated as before: 
U v (t + At/2) = U v (t) + (V-P) At/2 

3. V(t + At/2) = V(t) + Q- 1 ^^ + At/2) At/2 
4- r7 = r1(t) + I7I gi 7iy |At 

5. V(t + At) = V(t + At/2) + Q- l Tl v (t + At/2) At/2 
followed by two rescaling steps: 

. n{t + At) = ^r> 

# Pi ~ L(t+At)Pi 

6. V is evaluated using Eqn. |27j with r*j(t + At), L(t + At), and p"; then 
n v (t + At) = U v (t + At/2) + (V - P) At/2 

7. pi{t + At) = pi' + + At) At/2 . 

So far the algorithm has been developed for the case without friction and noise. For 
the case with friction and noise, we simply use the scheme given above, and introduce the 
following replacements: 

fiAt/2 -> fAt/2 - 7o— At/2 + JksTjoAtZi (28) 

(V - P) At/2 -+(V-P) At/2 - lv ^-At/2 + ^k B T lv Atz v . 

Here and denote uncorrelated random numbers with zero mean and unit variance; for 
simplicity we sample them from a uniform distribution via 

z = V12 (u - (29) 

where u is uniformly distributed on the unit interval. The momenta which occur in Eqn. 



28| are, for simplicity, taken as Pi(t) in step (1), ITy(t) in step (2), IIy(t + At/2) in step (6), 



and p'l in step (7). 
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B. Choice of Parameters 



We start from the observation that a molecular system is characterized by a typical 
molecular frequency c^o, which can be viewed as the "Einstein" frequency of oscillations of 
an atom in its "cage" |23|]. With use of the intermolecular potential v(r) one gets 

ul^^j d d rg(r)V 2 v{r), (30) 

with p being the particle number density, m the mass of the molecules, and g{f) the pair 
distribution function. Alternatively, one can define a molecular time scale by the time 
which a sound wave needs for traveling the nearest neighbor distance ||. However, both 
frequencies coincide by order of magnitude. This frequency governs the time step At which 
one has to choose in order to keep the MD algorithm stable; a typical rule of thumb says 
At = (l/50)(27r/wo). 

Similarly, the piston degree of freedom performs oscillations, if it is simulated with very 
weak friction in the NVT ensemble. Following Nose ||, we can estimate their frequency f2 
quite easily. Within a linearized approximation, the isothermal compressibility 

^-hw-vh^)-^ (31) 

controls the relation between pressure fluctuations SV = V — P and volume fluctuations 
SV = V - (V) via 

= w sv = -v^ 5v - (32) 



Therefore, one concludes from Eqn. [16] by trivial insertion 

d 2 _ 1 



SV = -——6V, (33) 
dt 2 QVk t ' v ; 

which is the equation of motion of a harmonic oscillator with frequency 

Obviously, the piston mass has to be chosen small enough such that the system can adjust 
its volume sufficiently fast. On the other hand, it cannot be chosen too small, since then Qq 



becomes too large, see Eqn. [34]. Clearly, one does not want f2 to exceed uo, since otherwise 
the simulation would need an unnecessarily small time step. The optimum piston mass is 
thus found from the resonance condition Q = uj |§, which yields a quite different value 
for Q (by a factor of L~ 2 ) than Andersen's original suggestion || - this original criterion 
has turned out to be not correct. The similar frequencies of the molecular oscillator and the 
box volume lead to a very quick energy transfer between them, resulting in a very efficient 
equilibration. However, one will often choose a substantially larger value for Q in order to 
separate the time scales, such that the molecular motion on short length and time scales is 
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largely unaffected by the piston motion. Regardless of the precise choice for Q, one should 
note that keeping Q constant implies a scaling of Q with the inverse system size. 

When the coupling to the heat bath with friction and noise is added, the question arises 
how to choose the damping parameters 70 and jy. Let us hence study a (deterministic) 
damped harmonic oscillator, 

mx + 72 + mWfli = 0. (35) 

Obviously, for small 7 the damping can practically be neglected. On the other hand, for 
7 « mouo, damping force and harmonic force are of the same order of magnitude (for a 
harmonic oscillator, we can estimate the velocity via x ~ ujqx). The exact calculation yields 

7 C = ImujQ (36) 

for the "aperiodic limit". This value quantifies the qualitative distinction between "weak" 
and "strong" damping, denoting the boundary between oscillatory behavior (7 < 7 C ) and 
pure relaxational dynamics (7 > j c ). Only in the weak damping case 7 <C 7 C , the fastest 
time scale (which governs the time step) is given by Uq, while for 7 ^> j c the damping term 
dominates, requiring a smaller time step. For this reason, 7 values beyond 7 C are clearly 
undesirable. Even worse, for 7 > ) c the relaxation contains also a very slow component, 
whose characteristic time is, for large 7, given by 7/(mwg). Therefore, 7 = 7 C is clearly the 
optimum value for fast equilibration. 

However, for the single-particle damping 70 we typically choose a value which is between 
one and two orders of magnitude smaller than 7 C . This is in accord with the philosophy of 
simulating the system rather close to its true Hamiltonian dynamics, such that at least on 
the local scales, both spatially and temporally, the dynamics can be considered as realistic. 
We hence use the coupling to the heat bath mainly for additional stabilization, deliberately 
keeping the molecular oscillations in the simulation. 

On the other hand, there is no analogous argument of "physical realism" for the box 
motion, which is intrinsically unphysical. One is therefore clearly led to the choice 7y = 7y c , 
such that the ringing is just avoided, while the box motions are still sufficiently fast. The 
same conclusion was obtained by Feller et al. |16| , where however no theoretical background 
was provided. 

Using this choice, one is left with only one time scale for the piston motion, given by Qq, 
and this is in turn adjusted according to the needs of the simulation: If one is only interested 
in statics, then "resonant" coupling is desirable (i. e. Qo ~ co> , and also 70 ~ 7 C ), at the 
expense of distorting the motions even on the molecular scale. If, on the other hand, it is 
desired to realistically simulate the molecular oscillations, one should enforce a separation 
of time scales by choosing both Q <C u and 70 <C 7 C . 

C. Numerical Test 

We study a system containing 100 particles interacting via a truncated Lennard- Jones 
potential whose attractive part is cut off: 
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U LJ (r) = 4e 



a\ 12 /a\ 6 1 
+ 



rj 4 



r < 2 1/6 (j. (37) 



We choose Lennard- Jones units where the parameters a and e as well as the particle mass 
m are set to unity. 

Figure [I] illustrates the problem of the "ringing" in this particular system, by displaying 
the autocorrelation function of the pressure fluctuations for various simulation parameters. 
The simulated temperature is k^T = 1.0, while the pressure was fixed at the rather low 
value P = 1.0. The mean volume for the 100 particles is V = 262.7, corresponding to 
a density p = 0.38. The compressibility at this state point is k,t — 0.3, such that for a 
box mass of Q = 0.1 one finds Qq ~ 0.36 for the ringing frequency or 2ir/Qo pa 18 for the 
oscillation time. The figure shows that the box indeed oscillates, and that the frequency of 
the oscillations has been estimated correctly. The left part of the figure is for pure undamped 
Andersen dynamics where the dependence on Q is shown. Choosing a value of Q = 0.0001 
leads to a very fast relaxation of the pressure autocorrelation function. The theoretical 
prediction for the best box mass for the molecular frequency of ujq pa 8.5 in this system 
leads to Q opt = 1/(VktuJq) ~ 0.00018. But as seen in Fig. Q, even for this value of Q the 
oscillations still remain on a very short time scale. Only a value of 7y = 0.001, close to 
7v c = 2QQ Q pa 0.002 suppresses the fluctuations efficiently and the autocorrelation function 
resembles the autocorrelation function in the NVE ensemble. Interestingly, it is also seen that 
for constant volume the pressure relaxation is considerably slower if the molecular damping 
is turned on. Conversely, this behavior is practically absent in the constant pressure case, 
where actually the "best" autocorrelation function was found for = 0.001 (as discussed), 
combined with some additional molecular friction 7 = 0.5. 

The statistical accuracy of an observable is given by the ratio between simulation time 
and the integrated autocorrelation time r, i. e. the value of the time integral over the 
normalized autocorrelation function [P4}| . From that perspective, a slow decay with many 



oscillations is actually not particularly harmful, since the integral value is rather small, 



due to cancellation. However, the result of Ref. |24| holds only in the asymptotic limit 
where the simulation time is substantially longer than the decay of the correlation function. 
Moreover, the numerical integration of an oscillatory function converges only slowly and is 
hence rather awkward. For these reasons, a simulation algorithm which avoids oscillations 
is clearly preferable. In order to illustrate this further, Table | lists r for the autocorrelation 
functions shown in Fig. |2|. The smallest r is actually found for the pure Andersen NPH 
simulation. However, turning on the damping increases r only by less than a factor of two, 
while the decay time decreases nearly by a factor of ten. 

The case of a box whose mass has, for reasons of separation of time scales, been chosen 
substantially larger than the optimum value, is illustrated in the upper right part of Fig. [I] 
(Q = 0.1, i. e. three orders of magnitude larger than Q opt ). The autocorrelation function 
decays most rapidly for j v = 0.1, as theoretically expected (jvc = 0.072). Compared to 
undamped dynamics at the same mass, this is a considerable improvement. Nevertheless, 
this decay is still substantially slower than what one can obtain if also the mass is chosen 
optimally (Fig. |2|) — this is simply the price which is being paid for achieving realistic 
molecular motion. For practical purposes, the decay obtained for Q — 0.1, 7y = 0.1 is quite 
acceptable. 
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VI. CONCLUSION 



We have discussed the algorithm of Andersen |Sf and the generalization to stochastic pis- 



ton motion by Feller et al. [16], generalizing it even further to also include stochastic motion 



of the particles. We gave a straightforward proof that the NPT ensemble is produced. The 
implementation by means of a symplectic algorithm is particularly stable and well-suited for 
MD problems. Another important point, the choice of the right simulation parameters, was 
studied both theoretically and numerically, and a guideline for their optimum values was 
given. We view this algorithm as a particularly useful realization of the constant-pressure 
ensemble. 
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FIGURES 
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FIG. 1. The left hand side shows the autocorrelation function of the instantaneous pressure 
fluctuations using the Andersen algorithm without stochastic forces. Various values of Q are used 
as indicated in the figure. On the right hand side the same function is displayed for stochastic 
dynamics, using the box masses Q = 0.1 and Q = 0.0001 and various damping constants "fv as 
indicated in the figure (molecular damping 70 = in all cases). 
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FIG. 2. Same as Fig. ||, but on shorter time scales. Even for the optimum value Q = 0.0001, 
box oscillations still remain, which do not occur in the pressure autocorrelation function in a 
constant volume simulation (right hand side). The use of a friction ~fy damps out the oscillations, 
and a combination of jy = 0.001 with 70 = 0.5 leads to the best agreement with the NVE result. 
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TABLES 



Q 


1 v 


/u 


r 


NVE 


NVE 


0.0 


0.042 


NVT 


NVT 


0.5 


0.127 


0.0001 


0.0 


0.0 


0.018 


0.0001 


10~ 5 


0.0 


0.023 


0.0001 


10~ 3 


0.0 


0.03 


0.0001 


10~ 3 


0.5 


0.031 



TABLE I. Integrated autocorrelation time r for the parameter combinations of Fig. |2|. 
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